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1. Introdnctioo 


When solving time-dependent partial differential equations numehc^y, it is often desirable to use 
implicit methods for reasons of numerical stability. This is particularly true for parabolic equations where 
tte time st^ restriction is severe for explicit methods. In more than one H>aoe dimesion, however, impli- 
cit methods lead to linear systems of equations with complicated structure which cannot be solved effi- 
ciently by direct methods. 

Consequently, a wide variety of dimensional splitting or fractional st^ methods have been used. In 
these methods, a single multidimensional implicit time step is rqtlaoed by a sequence of steps, each of 
which is implicit in only one coordinate direction. In that direction, the equations can be solved along one 
line of gric^xjints at a time, giving a banded syston of equations which can be easily solved. The reader is 
referred to [2], [4], [10], [13] for an introduction to many of the methods used in practice. 

One difficulty which has sometimes pl,<eued these methods and frequently caused confusion amor.g 
users is the proper ^xdfication of boundary conditions for the intermediate solutions which arise between 
the individual steps of this sequence. The purpose of this paper is to show how the proper boundary con- 
ditions can frequently by determined quite simply and logically by considering the "modified equation" (or 
"model equation") corresponding to each finite difference approximation in the sequence. 

We will illustrate this technique bj' considering several fractional step methods for the heat equation 
in two space dimensions: the locally one dimensional (LOD) method, the alternating-direction implicit 
(ADI) methods of Peaceman-Radiford[ll] and Douglas-Gunn[3], and the approximate factorization (AF) 
method of Beam and W’armingfl]. General boundary conditions involving normal derivatives are handled. 

We have chosen a simple problem and well-known fractional step methods in order to ilVvistrate the 
technique as dearly as possible. The same ideas used hert; can be applied to other methods and more com- 
plicated ec]uations. The extension to methods involving more than two fractional steps is illustrated by 
considering the LOD method in three space dimrisions in Section 6. 

Although here we only discuss dimensional splitting methods, the same approach can be used to 
determine the correct intermediate boundrjy conditions for any additive splitting method. Other contexts 
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witere these ideas have already been applied indude the flitting of hyperbolic systems into subproblems 
with di^>arate wave speeds{8], [9], die splitting of convection-diffusion equations into hyperbolic and para- 
bolic subproblems[8], and ^littings of the incompressible Navier-Stokes equatians[7]. 

For most of this paper we wiU consider the two-dimensional heat equation, 

u, = u„ - Lu in ft = [0,1] X [0,1] (1.1) 

with boundary conditions on dft. Along the boundary z = = 0 or 1), for example, we will consider 

either the Dirkfalet condition 

= g(y,0 (1.2a) 

or the more general Robbins condition 

a(y.0“(jCiO'.0 + ^(yr^K(*iO'.0 = Tf(y.O- (i-2b) 

A fracdonal step method for this problem typically has the form 

= P,U”j (1.3a) 

= P^ir^. (1.3b) 

Here the P’s and Q’s are spatial difference operators, with and involving differences in only the x- 

and y-directions re^)ectively. (/^ rqrresents the numerical aj^oxiniaticsn to the solution u(x,,yy,f„). For 

simpkaty, we take equal me^ widths h in the x and y directions, so x, = and yj = jh for 

ij = 0,1, ...,A^, where h = 1/N. The time stq) is denoted by k, so r„ = n*, n = 0,1,.... 

The problem is to determine a^qntjpriate boundary conditions for the "intermediate solution” 
which arises in (1.3). We typically require boundary data b = 0, N (or (1.3a). 'Ve may also require 
b — 0,N, for (1.3b) if P 2 involves y-differences. 

Typically, (1.3) gives a second order accurate approximation (T'*, while the individual stqs (1.3a) 
and (1.3b) are not second order accurate methods for the origiiul equation. Hence the intermediate solu- 
tion iT is iKm[^y»cal and does not coire^xmd to the true solution at any intermediate time. This makes 
it difficult to determine a priori the proper boundary conditions. 
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The approach we will t«kt> conststs of determining the modified equation for (1.3a), the differential 
equation which wodd be solved to second order accuracy if we iterated with the scheme 

The modified equation is derived by rqjladng the difierence operators by expansions of differential opera- 
tors. This expansion can be truncated at an appropriate point to give a difioential equation. 

We denote the modified equation by 

«* = L‘u (1.^) 

wlwre I* is a differential operator involving only ^tial derivatives. Similarly, we can derive the modified 

equation for (1.3b): 

«;* = £'*«“. (1.4b) 

By truncating die expansions in eadi modified equation ^^opriately, we will have 

£=£*-»- L'\ (1.5) 

an additive splitting of die ori^nal operator. 

There is an obvious flitting of L into one dimensional operators, namely 

r = d\, L" = a*. (1.6) 

This is precisely the splitting used by the LOD method discussed in Section 2. The odtsr methods, how- 
ever, use more complicated ^littings d the c^ierator L. 

Once the modified equation (1.4a) has been determined, we can view (1.3a) as taking a single time 
stqi on the equation (1.4a) with initial conditions 

= u{xc^,t„) V t ft. (1.7) 

Here we assume t/JJ = If the solution to (1.4a), (1.7) is denoted by u{x;y,i) for r t„, then 

Iflj = «*(x,ov.t„+*) + 0(*") (1.8a) 

since (1.3a) is second order accurate on (1.4a). One approach to specifying boundary conditions for (/* is 

to simply raqiloy a one-dded, second order accurate, finite difference approximation to (1.4a). Then we 

are solving the same equation at the boundary as in the interior, but with a one-sided scheme instead of 
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the (presumably centered) scfaone (1.3a). 

The problem widi this approach is that it does not make any use of the spedHed boundary conditions 
(1.2). We can often do better by determining the boundary behavior of u'(x^,t) in terms of the (given) 
boundary behavior of u. If we succeed in determining to 0(A^) for points along the boun- 

dary, this gives us the proper q)ecification of iT. 

Actually, it usually suffices Co have boundary data which ate one order of accuracy lower 
than Che interior scheme, i.e., we only need to insure that 

U“j = u(xiO>j,t^+k) -I- 0(*r) (1.8b) 

at boundary points to maintain second order accuracy. For hyperbolic equations this is a result of Gustafs- 

son[6]. For the heat equation it is a simple consequence of linearity and the maximum principle. How- 
ever, in practice it is found that such boundary conditions can lead to a large increase in the error constant. 
Although the results are second order accurate as the nesh is refined, the errors on any particular mesh 
may be an order of magnitude larger than necessary. This can make a significant difference in the effort 
required to solve a multidimensional problem to a specified accuracy, and so we recommend using boun- 
dary conditions with <2(k^) accuracy whenever possible. Some examples of the resulting increase in accu- 
racy will be seen in later sections. 

To determine the behavior of u" along the boundary, consider a typical point along one of the boun- 
daries X = jTj,. Expanding u"(xi„y,t„+k) about time t„ gives 

= u*(xiO-,r„) + ku‘(xj,o',0 + + • (1.9) 

~ u“ + kL'u' + • • • . 

We obtain Af* by differentiating (1.4a) with respect to t and rq)ladng time derivatives of u* on the right 
hand side by L‘u, so that A/* involves only tpatuU derivatives. For the linear equation considered here we 
have M - (L*)^, but the same process can be applied even in nonlinear problems. 
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It is important that the final expression in (1.9) involves only ^tial derivatives of u* at time t„. (We 
make the convention throughout this paper that if a function occurs with no arguments on the right hand 
side of such an expression, the point (x^o'>t„) is assumed.) According to the initial conditions (1.7), 
u* ■ II at time t„ and so L*u* = L'u and similarly for any spatial operator (N.B. thb is valid only at time 
r„t). So (1.9) becomes 

“*(jfi Ji'n+k) = « + kL*u + « + ■ • • . (1.10) 

This can frequently be manipulated to yield an expression in terms of the original boundary data for « by 
using the original equation (1.1) to replace L* and AT by "tangential" operators involving only r> and y- 
derivatives. The form of the final expression obtained dqjends on the particular case in question. Rather 
than continuing in such generality, we will demonstrate the process on particular examples in the following 
sections. 


2. The locaify one dimensional (LOD) method 

The LOD method for the heat equation (1.1) takes the form 


(/ - ^kD^U-j = (/ + 

(2.1a) 

- |kD')t/ry"' = a + 

(2.1b) 


Here we use 


D^U^j = - 2U,j + «,_i^), Dp,j = - 2Uij + «,j_i). 

Since (2.1a) is the second order accurate Crank-Nicolson method applied to 

«; = ( 2 . 2 ) 
we can take this as the modified equation for (2.1a). The second stq) (2.1b) is the same in the >• 

direction, and so the LOD method corresponds to the splitting (1.6). 

Note that (2.1a) can be {qjplied dX j = b = Q, N ia order to obtain the values iflj, needed in (2.1b), 
However, we stiM need to specify boundary conditions for V^j, = 0,N. Consider a typical point (^iO-) 
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and expand h*(x^o' .<»+*) about to obtain 

u(x^^,t„+k) = + • • • 

= u* + ihi* + — Jt^u* + • ■ ' 

= K + fa< + —k^u • 

"^xr 2 xxxx 

«4iidi corresponds to (1.10). We now use the original equation (1.1) to solve for 


(2.3) 


and 


Id = 1/ 2jd + £1 

■•rxxt ••ff **mT 


This allows us to rewrite the boundary data (2.3) in terms of t- and y-derivatives of u along the boundary 


X = X*; 


u'(xi.y,t„+k) = u + k(u, - u„) + yAr(u„ - 2«^^. + + 


= - *“>y Y*^(“my " Xy) + 

= [/ - kdf. + + - • • ] M(XiO',^„+*)- 


(2.4) 


This can also be obtained by integrating (1.4b) backwards in time from see Section 6 toi U.... d'. 

If the original problem specified Diricfalet boundary conditions, u(x^^,t) ~ g(y,t), (2.4) immediately 
provides the tq^propriate boundary conditions for (/*. Based on (1.8a), we can take 

= «(yy»^+*) - *f,y0y.^„+*) + •J*^«>yyy()'y.^+*)- (2 5) 

With the Robbins condition (2.2b), appropriate boundary conditions are still easy to derive if a and p are 
constant, so (2.2b) is 

««(jfAxV.O + 3«x(jf*o'.0 = y(y,0- (2-6) 

Differentiating (2.4) with respect to x and taking the linear combination au“ + Pu" gives the boundary 


condition 
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= >(y.^+*) - * 7 n(yA+*) + -|*^,TO(y.^n+*)- ( 2 . 7 ) 

If a and 3 are not constant then the situation is more complicated since, for exan^>le, au^y + fiu„y ^ 7^. 

Based on the second line of ( 2 . 4 ), we can still obtain a condition of the form 

au + p< = y(y,t„+k) - i(au^y + p«^^.) -f _ 2au^^. - 2pu,„y) , (2.8) 

but now the derivatives of u(x/,^,i„) on the right hand side must be replaced by finite difference approxi- 
mations based on C^. Unless this is done carefully, numerical inst^^ties can arise. This is discussed in 
more detail in Section 4 . 

To demonstrate the valicfity of the boundary conditions derived here, we presort the results of two 
numerical experiments with Diridilet conditions. We first consider the prohloc (1.1) with exact solution 

u(x^,t) = sin(x /2 + y)e~^ ( 2 . 9 ) 

and specify boundary conditions obtained by evaluating ( 2 . 9 ) along the boundaries, e.g., along x = 0 , 

^(y>0 = sin(y)e"^ (2.10) 

We compare tte results obtained with three different choices of boundary conditions Ulj. The first of 

these is simply 

ir,j = g(yj,t„+kn). (2.11) 

This has been recommended by some authors, but comparing this to ( 2 . 5 ) shows that it differs by 0 {k). 

We thus expect this boundary condition to destroy the second order accuracy. This is confirmed by the 
results of Table 2 . 1 . Here we have displayed the errors at time t = 0.75 measured in both the L.- and 
L2-norms. We have also computed die L2*rate of convergence, obtained by comparing the errors with dif- 
ferent values of k (we always take k/h fixed as the mesh is refined and h = VN). The is defined 

by \og{e-je{) / log(ityifci), where e■^ and are the L2-errors obtained with time steps ki and *2 respectively. 

The second boundary condition used is 

Kj = giyj,t„^k) - kg„,0>,t,+k), (2.12) 

which was derived by Dwoyer and Thames[ 4 ] using the "method of undetermined functions.” In this 

epproach, one inserts undetermined constants in place of the unknown boundary values and combines the 



8 


stages into a one-stq) method. Requiring that the resulting me±od be second order accurate yields the 
boundary condition (2.12). Since (2.12) corresponds to the first two terms in (2.5), according to our 
comments concerning equation (1.8b) this should be sufficient to restore second order accuracy. This is 
also confirmed by Table 2.1. Finally, we have used (2.S) itself. Again the method is second order accu- 
rate, but the errors are reduced by an order of magnitude over those seen with (2.12). 

In our second experiment, we investigate the effect of various boundary conditions on a steady state 
solution obtained with the LOD method. Presumably the choice of boundary conditions will have little 
effect on the rate of convergence to steady state, but may have a significant impact on the accuracy of the 
resulting solution. To investigate this latter effect, we have taken initial data equal to an exact steady state 
solution to (1.1), namely, 

“(jf«y) = cosh(y~-5) sinx cosh(x-.5) siny. (2.13) 

Again, Diridilet boundary conditions were obtained by evaluating (2.13) along fiie boundaries. Within a 

few iterations, the LOD method converges to a numerical steady state. Table 2.2 shows the error in the 

numerical solution for each of the boundary conditions considered previously. Again, (2.11) gives only 

first order accuracy while (2.12) restores second order accuracy. In this case, the use of (2.5) improves 

the error constant by a factor of at least 20 over (2.12). 


3. The Pcaceman-Rachford and Dooglas-Gunn ADI Methods 

The ADI method introduced by Peaceman and Rachford[ll] has the form 


(/ - ^kDj)Ul = 

(3.1a) 

(/ - ^kDj)U”;^ = (/-»■ ^kD^U^j • 

(3.1b) 


Here the splitting is no longer simply (1.6). The ADI method is usually viewed in the following way: 
Each stq) (3.1a) and (3.1b) is a first order accurate method for the original equation (1.1) on a time stq) 
of 'ength k/2 which combine to give second order accuracy over a stq) of length k. Because (3.1a) is con- 
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slstent with the original equation, it is ten|)tiiig to specify 

iKj = (3.2) 

for example, in the case of Diridilet conditions. However, it has long been known that this causes a loss 

of accuracy since this boundary condition does not contain the 0(lr) error present in the interior which is 

required to cancel out the 0(h^ error in the second step. Fairweather and Mitchell[5] determined the 

correct boundary conditions by the method of undetermined functions, obtaining 

vlj = jQ - !(/ + ^kDj)gf (3.3) 

where gj - g(yj,t„). The same result can be obtained by the approach advocated here. 

Rather than viewing (3.1a) as a first order scheme for (1.1) with time step k/2, we derive the modi- 
fied equation for which (3.1a) is a second order {qjproximation with time step k. We have 

(/ - = (/ + ^kD^)u(x,y,t„). 

Expanding u"(x,y,t„+k) about «"(jt,y,r„) and using, for example, 

= ay + 0(k^) 


we obtain 

u + k(u‘ - y«;,) + Y*^(«“ - «"„) + 0(k^) = «“ + jku'^ + 0(k^) 

or, solving for 

+ «n) “ + 0(fc2), (3.4) 

We can obtain expressions for u’ and by differentiating this. Plugging the resulting expressions back 
into (3.4) gives 

= y(«« + “>>) + - «m>) (3.5) 

plus terms which are 0(Jr). Hence (3.1a) gives a second order accurate approximation to (3.5) on a time 
step of length k and we can take 
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L" = i(a2 + dp + jk(di - dp. (3.6a) 

Since (3.1b) is the same as (3.1a) but with x and y interchanged, deriving the modiHed equation for 
(3.1b) gives 

+ 5y) + |*(«v - ^t)- (3.6b) 

Note that (1.5) is satisfied. 

Using the equation (3.5), we can proceed as before to determine the boundary data. Since 

L" = + |-fc(dj - dp, hr = (L‘)2 = + o(k), 

we now obtain, by (1.10), 

= 14 + ^kLu + |-*(i“ - ^L)u + ^k^Ou + • • • (3.7) 

= - «m>) + • • • • 

Moreover, we can manipulate (1.1) to obtain 

U “ M M ” 2u 

(t *»■ 

so that this can be reexpressed in terms of t- and y- derivatives along the boundary. Hence, with Dirichlet 
boundaiy conditions, we can take 

Hi] = g(y],tn^fH2) + -|lr(g„(>'y,r„)-2g,,,.(yy,f„)). (3.8) 

This gives the required 6>(lr) correction to (3.2). Instead of using derivatives of the function g, one could 
approximate (3.8) by any finite difference expression correct to 0(Jlr^). It is easy to verify that 
Fairweather and Mitdiell’s boundary condition (3.3) is one sudi approximation. 

The approach taken here easily extends to more general boundary conditions. For example, 
corresponding to the boundary conditions (2.6) with a and p constant, we obtain 

a«“(x*,y,f„+*) + 3i4;(x*,y ,/„+*) = y(y,i‘„+k/2) + ■r*^('y„(y,0 “ 2y^.,.(y,f„)). 

© 
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Based on the formulas used above, it is also easy to determine the intermediate boundary conditions 
for the ADI method proposed by Douglas and Gunn[3]: 

(7 - kD^Ul = (/ + kD^U^j (3.9a) 

(/ - kDj)U”;^ - f/Ty + (3.9b) 

Except for the factor 1/2, (3.9a) is identical to (3.1a) and the modified operator is easily determined to be 

L‘ = (aj + a?) + i*(a? - d% (3.io) 

The intermediate boundary conditions in the Dirichlet case are found to be 

Vlj = ,.'„+*) + yfc2(«„(yy,fj - 2g^y(yj,t„)). (3.11) 

Although it is not necessary to con^te the operator L'*, it is interesting to do so since the equation 
(3.9b) has a feature not seen before; it involves U” as well as iT and Because of this it appears 

senseless to discuss the modified equation for (3.9b) since we cannot apply it in isolation. However, if we 
multiply (3.9b) by (/ + kD;) and use (3.9a), we can eliminate U” and obtain 

(/ + kDj)(l - kD;)U”^' - [(/ - kD^) + k(I + fc9;)D;] (T. (3.12) 

Note that this is a bit of a twist on the usual analysis in which one eliminates iT from (3.9) to obtain a 

second order accurate method for (1.1). We have eliminated V" to obtain a second order accurate method 

for the (as yet unknown) modified equation (1.4b). In this equation we use initial data 

«**(r„+*)-«'(r„+k) = fr 

and view C/"'*'* as an 0(*^) approximation to «““(f„+2Jt). Replacing by «““(f„+Jfc) + 

ku"’'(t„+k) + • • • ir (3.12) and proceeding as usual gives 

L" = |k(aj - di) 

after a tedious calculation Of course the same result can be obtained mudi more easily by using (3.10) 
and (1.5). However, this sort of technique is sometimes necessary in dealing with multi-s^ep methods 
involving more than one intermediate solution. 
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4. Hw Ai^wdmate Fttori a thm (AF) Method 

Ai^oxiinate factorization methods are typically written in a form different from (1.3), namely 

= Pir^ (4.1a) 

= K (41b) 

where 8" = - U" so that the new approximation is obtained by setting 


ur;^ = u?j + h:j ( 4 . 2 ) 

after solving (4.1). Qearly (4.1) is a two-stqp procedure for solving 

Q.QyV?;' = (^ + 

and Q^Q^. is an approximate factorization of some two-dimensional spatial operator. 

For the AF method, we need only specify boundary conditions for 6‘ along x = x^. If Diridilet con- 
ditions are imposed, then the proper boundary conditions for hlj in (4.1a) can be easily determined using 
(4.1b): 


Kj - QM, 

= Gv(*(yy.'n+i) “ «(yy.^))- 


(4.3) 


For more general boundary conditions, however, it is useful to pursue the modified equation approach. 
Beam and Warming[l] consider AF methods for a more general problem «vith mked derivatives: 


«C = + *«T> + (4-4) 

They show that stable second order accurate methods can be obtained while handlLig the mixed derivative 

term explicitly. One such method is 


(/ - a,*oD2)8- = ^ 

(/ - wkaDf)h”j = S,}. 


k(aD'; * cDf)(I + (C-0+|)V) kbDJ)^.(I + (?+y)V) -t- 




Here u, and 6 are parameters satisfying 


(4.5a) 

(4.5b) 


0 ) = 0/(1 0 
and 
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VU" = IT - = 8"~\ 

- t/(-w.: + t/.-v-:)- 

For the heat equation a = c = 1, 6 = 0, and the sin^>leat such mediod (corresponding to Crank-Nkolson) 

is obtained by taking ( = 0, o> = 6 = 1/2. Then the operator P on the li^t hand side of (4.5a) reduces 

to 

P = k{D^ + D^. 

For b * 0, other parameter choices give more stable methods, e.g., 6 = 1, ( = 1/2, o> - 2/3 (see [1]). 
We have also found that with Robbins boundary conditions, the latter choice of parameters gives a more 
stable method even for the heat equation. 

In this section, we will consider the more general method (4.S) since the additional complicatims 
only affect the operator P in (4.Sa). We will see diat this has essentiaUy no effect on the derivation of 
boundary conditions. 

Thus the method has the form (4.1) with 

(?, = (/- otkaD^) (4.6a) 

Q, = (/ - u)*aD;) (4.6b) 

and P given by (4.5a) (P could be even more complicated - e.g., for the incoic^essible Navier-Stokes 

equations we could obtain a method of this form in which P contains nonlinear terms [7]). In order to 

derive the modified equation for (4.1a) without having to expand all the terms in P, we instead base our 

expansion on (4.1b): 

QyiU?r^ - U?j) = l/Ty - U?j- (4.7) 

Assuming as usual that U“j = a" (x,,ypt„+k) with irutial conditions (1.6), we expand (4.7) to obtain 

(/ - <akcd^{ku + + ••) = *“,'+ + •■• 

all evaluated at time t„. Using u, = Im = Lu at time t„, solving for u* gives 
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«* = Lu" + ^k(L' — Tuacd^L - d^u'+... ^ (4.8) 

Differentiatiiig this with reflect to t gives ul — L^u" + 0(k). Inserting this in the final term of (4.8) 
>ields 

u‘ = Lu" — iakcd}Lu (4.9) 

plus terms wiiich are (7(it^. Thus (4.Sa) gives a second order accurate approximation to (4.9) if we set 

iT = U”+t". We thus have 

r = (/ - utkcdf)L. (4.10) 

We first consider the case of Diridilet boundary oonditioiis. Then (1.10) gives 

= « + *(/- (tikcdf)Lu + ^krL-u + 0(k^) 

since (L')~ - Lr + 0(k). Ufing the origina] equation u, = Lu, we can simplify this to obtain 

= u(xi„y,t„+k) - <akrcu,,y(Xi„y,t„) + 0(1^). f4.11) 

We obtain the bouiKlary condition for 6' by subtractLig u(x^^v,r„) - g(y,t„) from both sides, and dropping 

the 0(fc’) terms, 

Kj = - 8(ypfn)) - ('*-12) 

Notice that using sinqjly 

K) = giypfn^:) - g(yjfO ('^-i^) 

should give second order accuracy, but again including the final term in (4.12) can lead to a significant 

improvement in the error constant. It is interesting to compare (4.12) to the bouiKlary conditions (4.3) 

derived directly from (4.1). Using (4.6b), it is easy to check that they agree to 0(k^). 

Table 4.1 shows some numerical results with Diricfalet boundary conditions for the same problon 
used in Section 2. We see that (4.10) does give second order accuracy, but that the results are improved 
by using (4.9). Note that for a steady state problon with time-independent boundary conditions, - 0 
so that (4.10) and (4.9) agree. 
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For the Robbins boundary conditions (1.2b), we can differentiate (4.11) with rennet to x and take 
ths appropriate linear combination to obtain 

+ ^""V(x*0',^,+*) = r*' - (4.14) 

where we use the shorthand a""^- - a(y,t,»i), for example. Subtracting a"*'u(xj,jf,t„) + 3"" «,(jr*j',rJ 

from (4.14) gives boundary conditions for 8*; 

a-' 6' + 3-*-S; = y-*' - (a-^'ii + - «*=c(a"^ «„. + P"**^) (4.15) 

= (y"^- - y") - ((a-*'^ - a-)u + (P-- - p-)uj - «*=c(a"*'u„, + 3"^ -«„>,)• 
Here we are abusing notation slightly and lettiqg 8* also rqrasent the function u"(xj,t) - u(x^,t„) 

evaluated at r = Note diat if a and p are constant, (4.15) becomes simply 

a8* + P8: = (y"*- - y") - <-ilrcy-,. (4.16) 

When a and p are not constant, it is necessary to discretize the derivatives of u(xj,^,t„) occurring on 
the right hand side of (4.15). This can easily be done to the required accuracy, but Li practice it appears 
that great care must be taken in order to avoid numerical instabilities caused by high order differences in 
the boundary conditions. One possibility is to transfer as many derivatives as possible onto the functions 
a, p and y, as we now discuss. First note that if we assume the nond^eneracy condition P(y,r) # 0 for 
all y, r, we can divide (1.2b) by P(y,r). Hence, by modifying a and y we can assume without loss of gen- 
erality that P “ 1. 

With this assumption, coitqMting gives 

yn> = + “r) 

= -i- + o„,.H -t- ajyU, + 2a„u,. + 2o,«„ + o,«„. 

Solving for au^^. + gives ar. expression which can be used on the ri^t hand side of (4.15) (note that 

= o"u^y + 0(k)). Then (4.15) becomes 

a"*'8‘ + 8; = (y"^‘ - y") - (a"*' - o")« (4.17) 

- - o”yU, - 2a>^. - 2a>,j. - a>^y]. 

This can be discretized in a strai^tforward manner and appears to give nuch more stable boundary 
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conditions than directly discretizing (4.15). For oomfdeteness we present the details of our AF implemen- 
tation in Section 5. 

As a mimerical eq)eriment, we compare (4.17) to the more naive conditions 

= (Y- - y”) - - a")u (4.18) 

which do not oontiun the 0(k^ correction. We again use the heat equation (1.1) with solution 

«(x^,r) = ^x-^y-l) (4-19) 

For simplicity, we ^lecified Diikfalet boundary conditions along three boundaries and the Robbins condi- 
tion 

a(y,r)u(ljf,t) + u,(lj».f) = y(y,r) 
only along the boundary jc = 1. We take 

o(y,r) = costy) e', y(y,t) = Ysin(2y) e"' + cos(y) 

Table 4.2 gives a comparison of result obtained at r = 0.3 with three different combinatiom of boundary 
conditions: 

a) (4.13) at = 0 and (4.18) at x = 1, 

b) (4.12) at X = 0 and (4.18) at x = 1, (4.20) 

c) (4.12) at X = 0 and (4.17) at x = 1. 

As expected, all three ocnnbinations give second order accuracy, but including the correct 0(t^ terms 
improves the error constant. The combination (4.20c) gives errors about 10 times smaller than (4.20a). 
In this experiment we have used the parameters 6 = 1, | = 1/2, « = 2/3 for improved stability. 

5. I mp temen tattoo of ttie AF Method 

We include a detailed discussion of our implementation of (4.5) since in performing the experiments 
presented above it was found that both the accuracy and statnlity of the scheme were greatly aff cted by 
the manner in which the boundary conditions were imposed. 
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For each fixed j - 1, 2, A’-l, (4.5a) gives rise to a tridiagonal system of equations of the form 


~ (1 + 2uira)h'j ~ (i)ra8,>;^ = p^j (5.1) 

for t = 1, 2, A^-1, where r = k/hr and is given by the right hand side of (4.5a). Note that 

VUfj - 8,}~- which can be saved from the previous time step. Since in general this is a three level 
scheme, two time levels are required as initial data. In our numerical experiments we used the exact solu- 
tion at tine 0 and k, but in i ractice a two-level .scheme must be used for the first step. 

When Dirichlet conditions are imposed, the system (5.1) is ocxnpleted by specifying and Syy as 
disaissed in Section 4. For the Robbins condition we discretize the boundary conditions obtained for 8* in 
Section 4. For example, when a and ^ are constant we disaetize (4.16) as 

<^Kj - Pbj (5.2) 

for 8 = 0 or N, with rq)resenting the right hand side of (4.16). This can be combined with equation 

(5.1) for i = i to eliminate 81, (when h = 0) or 8y^.^ (when b = N). For example, at the right boundary 
(fc = N) we obtain 

^ f-(P.N7 - «8;y) (5.3) 

from (5.2). Inserting this in (5.1) gives the equation 


— 2wra8v_.y [1 + 2tara(l + Ao/p)] 8yy = py/ 

with a similar equation obtained for i = G. We thus have a tridiagonal system of N 1 equations for 8,y, 
i = 0, 1, ...,N. 


Note that the expiession for p,y (the right hand side of (4.5a)) involves second differences of Ufj and 
6,y“- in the jc-direction and henoe evaluating pyy requires values f/y-i,/ These can be obtained 

in the same way we derived (5.3) by discretizing (1.2b). We obtain 

= Uy-ij + y ( t; - “t^v;) (5.4) 


and 
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Actually, in practice we have found that the use of (5.4) leads to a considerable loss of accuracy 
near the boundary vMdi can be av*oided by using a hi^r order approxination to the derivtivf^ If 
(1.2b) is instead discretized by 

+ mj - m-'.j + v”s..j) = >; 

we obtain 

- V^s-ij) + j(*y; - at^y) (5.5) 

in place of (5.4). This did not seem to adversely affect the numerical stability but in^oved the errors 
considerably. 

When a dqiends of >■ and t (and ^ * 1 as discussed in Section 4), the boundary condition (4.17) is 
discretized to again give (5.3) with a replaced by and 

Psj = (nV' - T") - (“"" - - <yU”si - - V”sj-^yh 

- mj + 


6. Mnltl-stiep methods 


In many situations it is necessary to use fractional stq) methods involving more than two steps. Con- 
sequently, additional intermediate solutions arise for which boundary conditicns must also be specified. 


For exanq)le, the LOD and />F methods are easily extended to three space dimensions by adding a 
third stq} to the process. We consider the LOD method, whidi becomes 


(/ - ^kD^ir = (/ -H 

(6.'ta) 

(7 - ^kDf)U“ = (I + 

(6.1b) 

(7 - ^kDt)U"^^ = (7 + ^kD^V". 

(6.1c) 


For darity we have dropped the (three) subscripts on these variables and will also siq)press the spatial 
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argument of funcdom where convenient. The method (6.1) corresponds to the splitting 

L = L* + L** + L“* 

with C and I“* given by (1.6) and V" = d]. The boundary conditions for u along x = are again 
given by (2.3), but now we have 

“« = (i- - f-" - = «, - - “a 

= (L - L" - L-fu. 

The boundary conditions for tT along x = Xf, become 

Ulij=[I - k(dj + 3^ + »„+*:) . (6.2) 

Note that this involves only tangential derivatives along the boundary x = x^. 

In the second stq), we need boundary data for (/" along y - yi,. Recall that we view iT’ as an 
a^jproximation to u“(t„+2k), where «'* satisHes (1.4b) with initial coiKhtions u““ ■ m" at r = t„+k. 
Hence we can expand 

u\t„-^2k) = [/ + kL“ ^ ^k^(L"y + ■■■] u{t„+k). (6.3) 

Using the previously determined expansion (1.10) iii (6.3) gives boundary data for lT‘ in terms of u(t„). 

Actually, for three-step methods such as this, there is a miuh easier approach. By (6.1c), we can 
view {T* as an ^Jproximation to the function obtained by solving (1.4b) backwards in time from 
U"'*'’ = u(t„+.) With time step -* we thus obtain 

ir- = [i- kL- + |*2(l**') 2 - ...] «(f„,i) 

so that the boundary conditions are 

V‘^j = u~ ku„ + 

all evaluated at (xi^^,Zj,t„^y 

Finally, we note that even in two space dimensions it may be necessary to use a three-step method if 
the operator L is split as the sum of two noncxmunuting operators, as is typcally the case in variable coef- 
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fident or nonlinear problems. Then second order accuracy is retained by using a three-step Strang split- 
ting[12]. 
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T'«ble2.1 


BC 

N 

a 

L«-crror 

Li-enot 

L2-ratc 


11 

.15 

8.14E-3 

3.69E-3 


(2.11) 

21 

.075 

4.50E-3 

1.97E-3 

.91 


31 

.05 

3.17E-3 

1.35E-3 

.93 


11 

.15 

1.54E-3 

7.05E4 


(2.12) 

21 

.075 

4.38E4 

1.94E-4 

1.86 


31 

.05 

2.04E-4 

8.93E-5 

1.91 


11 

.15 

1.07E-4 

5.49E-5 


(2.5) 

21 

.075 

1.51E-5 

8.66E-6 

2.66 


31 

.05 

4.73E-6 

3.04E-6 

2.58 


Table 2.2 


BC 

N 

a 

Z.^-error 

Li-enor 

1 , 2 -rate 


11 

.15 

1.55E-2 

5.28E-3 


(2.11) 

21 

.075 

9.87E-3 

3.02E-3 

.81 


31 

.05 

7.23E-3 

2.11E-3 

.88 


11 

.15 

2.64E-3 

1.06E-3 


(2.12) 

21 

.075 

7.94E-4 

2.93E-4 

1.86 


31 

.05 

3.76E-4 

1.34E4 

1.93 


11 

.15 

1.02E-4 

4.45E-5 


(2.5) 

21 

.075 

1.59E-5 

8.59E-6 

2.37 


31 

.05 

6.09E-6 

3.41Et6 

2.28 
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Tabic 4.1 


BC 

N 

a 

L,-crror 

L2*aTor 

L2*rate 


11 

.15 

2.12E-3 

9.55F--4 


(4.10) 

21 

.075 

5.75E-4 

2.50E-4 

1.93 


31 

.05 

2.63E-4 

1.13E4 

1.96 


11 

.15 

8.65E-4 



(4.9) 

21 

.075 

2.57E4 


1.80 


31 

.05 

1.20E-4 

3.64E-5 

1.88 


Table 4.2 


BC 

N 

D 

L,-error 

£. 2 -error 

t2-rate 


11 

.06 

1.35E-3 

4.38E-4 


(4.20a) 

21 

.03 

4.11E-4 

1.15E4 

1.93 


31 

.02 

1.94E-3 

5.20E-5 

1.96 


11 

.06 

6.55E-4 

2.41E4 

■ 

(4.20b) 


.03 

1.64E-4 

5.57E-5 



B 

.02 

7.32E-5 

2.41E-5 



11 

.06 

1.40E-4 

4.96E-5 


(4.20c) 

21 

.03 

3.46E-4 

9.80E-6 

2.34 


31 

.02 

1.54Er5 

4.00E-6 

2.21 
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